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Abstract 

A two space dimensional active nonlinear nonlocal cochlear model is formulated in 
the time domain to capture nonlinear hearing effects such as compression, multi- 
tone suppression and difference tones. The micromechanics of the basilar membrane 
(BM) are incorporated to model active cochlear properties. An active gain parame- 
ter is constructed in the form of a nonlinear nonlocal functional of BM displacement. 
The model is discretized with a boundary integral method and numerically solved 
using an iterative second order accurate finite difference scheme. A block matrix 
structure of the discrete system is exploited to simplify the numerics with no loss 
of accuracy. Model responses to multiple frequency stimuli are shown in agreement 
with hearing experiments. A nonlinear spectrum is computed from the model, and 
compared with FFT spectrum for noisy tonal inputs. The discretized model is effi- 
cient and accurate, and can serve as a useful auditory signal processing tool. 

Key words: Auditory signal processing, cochlea, nonlinear filtering, basilar 
membrane, time domain 
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1 Introduction 



Auditory signal processing based on phenomenological models of human per- 
ception has helped to advance the modern technology of audio compression 
[10]. It is of interest therefore to develop a systematic mathematical frame- 
work for sound signal processing based on models of the ear. The biomechanics 
of the inner ear (cochlea) lend itself well to mathematical formulation ([1,2] 
among others). Such models can recover main aspects of the physiological 
data [13,11] for simple acoustic inputs (e.g. single frequency tones). In this 
paper, we study a nonlinear nonlocal model and associated numerical method 
for processing complex signals (clicks and noise) in the time domain. We also 
obtain a new spectrum of sound signals with nonlinear hearing characteristics 
which can be of potential interest for applications such as speech recognition. 

Linear frequency domain cochlear models have been around for a long time 
and studied extensively [8,9]. The cochlea, however, is known to have nonlinear 
characteristics, such as compression, two-tone suppression and combination 
tones, which are all essential to capture interactions of multi-tone complexes 
[4,5,17]. In this nonlinear regime, it is more expedient to work in the time 
domain to resolve complex nonlinear frequency responses with sufficient ac- 
curacy. The nonlinearity in our model resides in the outer hair cells (OHC's), 
which act as an amplifier to boost basilar membrane (BM) responses to low- 
level stimuli, so called active gain. It has been shown [3] that this type of 
nonlinearity is also nonlocal in nature, encouraging near neighbors on the BM 
to interact. 

One space dimensional transmission line models with nonlocal nonlinearities 
have been studied previously for auditory signal processing [5,16,17,15]. Higher 
dimensional models give sharper tuning curves and higher frequency selectiv- 
ity. In section 2, we begin with a two space dimensional (2-D) macromechani- 
cal partial differential equation (PDE) model. We couple the 2-D model with 
the BM micromechanics of the active linear system in [9]. We then make the 
gain parameter nonlinear and nonlocal to complete the model setup, and do 
analysis to simplify the model. 

In section 3, we discretize the system and formulate a second order accurate 
finite difference scheme so as to combine efficiency and accuracy. The matrix 
we need to invert at each time step has a time-independent part (passive) and 
a time-dependent part (active). In order to speed up computations, we split 
the matrix into the passive and active parts and devise an iterative scheme. 
We only need to invert the passive part once, thereby significantly speeding 
up computations. The structure of the system also allows us to reduce the 
complexity of the problem by a factor of two, giving even more computa- 
tional efficiency. A proof of convergence of the iterative scheme is given in the 
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Fig. 1. The figure on the left is a schematic of the cochlea, while the figure on 
the right represents the upper chamber with the macromechanical equations and 
boundary conditions. 



Appendix. 



In section 4, we discuss numerical results and show that our model successfully 
reproduces the nonlinear effects such as compression, multi-tone suppression, 
and combination difference tones. We demonstrate such effects by inputing 
various signals into the model, such as pure tones, clicks and noise. A nonlinear 
spectrum is computed from the model and compared with FFT spectrum for 
the acoustic input of a single tone plus Gaussian white noise. The conclusions 
are in section 5. 



2 Model Setup 



2. 1 Macro-mechanics 



The cochlea consists of an upper and lower fluid filled chamber, the scala 
vestibuli and scala tympani, with a shared elastic boundary called the basilar 
membrane (BM) (see Figure f). The BM acts like a Fourier transform with 
each location on the BM tuned to resonate at a particular frequency, ranging 
from high frequency at the basal end to low frequency at the apical end. The 
acoustic wave enters the ear canal, where it vibrates the eardrum and then 
is filtered through the middle ear, transducing the wave from air to fluid in 
the cochlea via the stapes footplate. A traveling wave of fluid moves from the 
base to the apex, creating a skew-symmetric motion of the BM. The pressure 
difference drives the BM, which resonates according to the frequency content 
of the passing wave. 

We start with simplification of the upper cochlear chamber into a two dimen- 
sional rectangle Q = [0, L\ x [0, H] (see Figure 1). Due to the symmetry, we can 
ignore the lower chamber. The bottom boundary (z = 0) is the BM, while the 
left boundary (x = 0) is the stapes footplate. The macromechanical equations 



3 



are 



Ap(x,z,t) = ^ + ^ = 0, xe[0,L}, ze[0,H], te[0,oo) 
p x (0,z,t)=T M p c (t), p(L,z,t) = (2.1) 
p z (x, 0, t) = 2pu tt , p z (x, H, t) = 

where p(x,z,t) is the pressure difference across the BM, u(x, t) denotes BM 
displacement, and p is fluid density. 

At the stapes footplate (0,2;), p c {t) is pressure at the eardrum while Tm is a 
bounded linear operator on the space of bounded continuous functions that 
incorporates the middle ear filtering characteristics. In the frequency domain, 
for each input e luJt , T M (u) = 2piuj/Z m , where Z m is the impedance of the 
middle ear. The middle ear amplification function is given by om = |2~m|. In 
our case, based on Guinan and Peake [7], 
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a M (f) = 1.815/ 2 ((1 - + (2C m /// m ) 2 )- 1/2 , (2.2) 

J m 

where / m = 4 kHz is the middle ear characteristic frequency and Cm = 0.7 is 
the middle ear damping ratio. Thus, for p c (t) = A exp{27ri/} + c.c, where c.c. is 
complex conjugate, we have T M p e (t) = B exp{27ri/} + c.c, where B = ay[(f)A. 

For more complex stimuli, it is useful to model the middle ear in the time 
domain as a one-degree of freedom spring-mass system. The equivalent time 
domain formulation of the steady state middle ear is given by 

p e (t) = m m s(t) + c m s(t) + k m s(t) ^ 
s (0) = s(0) = 

where s(t) is stapes displacement and m m , c m and k m are the mass, damping 
and stiffness of the middle ear. The stapes boundary condition in (2.1) is 
replaced by 

Px (0,z,t) = 2ps(t) (2.4) 

One of the interesting effects of using the time domain middle ear model is 
that it reduces the dispersive instability in the cochlea (see [14]). It appears 
that the steady state middle ear model ignores important transient effects and 
phase shifts that help to reduce the shock to the cochlea. 

At the helicotrema (L,z), we have used the Dirichlet boundary condition 
p(L,z,t) = 0. In [9], they used an absorbing boundary condition p x {L) = 
cpt(L), where c is a positive constant [9]. Other models use the Neumann con- 
dition p x (L) = 0. It has been stated that the frequency domain solutions are 
minimally affected by which boundary condition is chosen [8], and thus we 
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Fig. 2. Cross section micromechanics of the 
cochlea. The mass mi represents a cross sec- 
tion of the BM, while mass m 2 is a cross sec- 
tion of the TM. (Reconstructed from Figure 3 
in [9]) 



have chosen the simpler Dirichlet condition. However, interesting results on 
choosing the best initial conditions to minimize transient effects (dispersive in- 
stability) has been shown in [14] using the Neumann condition. To summarize, 
the macromechanics consist of equations (2.1)-(2.4). 



2.2 Micromechanics 



This is a resonant tectorial membrane model based on [9]. The BM and TM 
(tectorial membrane) are modeled as two lumped masses coupled by a spring 
and damper, with each mass connected to a wall by a spring and damper. 
(See Figure 2). A classical approximation is to have no longitudinal coupling 
except that which occurs through the fluid. Denoting t) = (u(x, t),v(x, t)) 
as BM and TM displacement, respectively, the equations of motion for the 
passive case at each point along the cochlea are given by 



M v i + C v i + K v i = F 



where 



mi 


) Cp — 


m 2 





C\ + c 3 -c 3 
-c 3 c 2 + c 3 



ki + k 3 -k 3 
-k 3 k 2 + k 3 



(2.5) 



(2.6) 



and forcing function 



(2.7) 



p(x,0,t) 


The parameters rrii, Cj, and fcj are functions of x. The initial conditions are 
given by 

S(x,0)=£(x,0) = (2.8) 
To make the model active, a self-excited vibrational force acting on the BM 
is added to (2.5): 

M p £ + C P i + K p £ = F + F a 
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where 



7[c 4 (m — v) + k 4 (u - v)\ 




The difference u — v represents OHC displacement. The parameter 7 e [0, 1] 
is the active gain control. In [9], this is a constant, but in our case will be a 
nonlinear nonlocal functional of BM displacement and BM location. Bringing 
F a to the left, we have 



where 



M p £ + (C p - 7 C a )e + (K p - 7 K a )e = F 



a, = 



C4 — C4 






— &4 


















(2.9) 



(2.10) 



Thus, the micromechanics consist of equations (2.6)-(2.10). 



^.5 Nonlinear Nonlocal Active Gain 



A compressive nonlinearity in the model is necessary to capture effects such 
as two-tone suppression and combination tones. Also, to allow for smoother 
BM profiles, we make the active gain nonlocal. Thus we have 

u(x,t) = f L e-^ 2 / x u 2 (s,t) s 

V Xrr Jo 

and gain 



where 6, X are constants 



7(M) = IT^ 



2.4 Semi-discrete Formulation 



Solving the pressure Laplace equation on the rectangle using separation of 
variables, we arrive at 

00 

p(x,0,t) = T MPe (t)(x -L) + J2 A n cosP n x (2.11) 

71=1 

where 




(2.12) 



6 



(n- |)tt 



L 



Substituting (2.11) into (2.7) and then discretizing (2.9) in space into N grid 
points, we have 

Mg u + C(t)g t + K(t)Z = b(t) (2.13) 

where 

Mi + «M f 
M 2 



M 



c(0 = c p - f(t)c a 



AT(t) = K p - f(t)K a 



Cl + c 3 -r(t)c 4 -(c 3 -r(*)c 4 ) 
+ x 3 - r(t)tf 4 -(k 3 - r(t)K 4 ) 



b(t) 



-K 3 

T M p e (t)(x - L) 




K 2 + K 3 



K 



k=l 



coth f} k H 



a 



cos(f3 k Xi) cos(f3 k Xj) 
ApH 



w ; 



N — 1 



C p , K p , C a and K a are now block diagonal, where Ki = diagjfcj} and Cj = 
diag{q}. Also, Mj = diagjm,}, T(t) = diag{ 7i (t)} and f(t) = diag{r(t), 0}. 
The numbers are numerical integration weights in the discretization of 
(2.12) and are chosen based on the desired degree of accuracy. Note that we 
can write Mf = M^W, where W = diag(wj) and Mf is symmetric and positive 
definite. The result of separation of variables produced the matrix Mf, which 
is essentially the mass of fluid on the BM and dynamically couples the system. 



3 Numerics 



In formulating a numerical method, we note that the matrices in (2.13) can be 
split into a time-independent passive part and a time-dependent active part. 
In splitting in this way, we are able to formulate an iterative scheme where 
we only need to do one matrix inversion on the passive part for the entire 
simulation. Thus, using second order approximations of the first and second 
derivates in (2.13), we arrive at 

(L p - L'2)i n+1 = B n = L p 1 B n + L p l L n J n+1 > k (3.14) 
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where superscript n denotes discrete time, k denotes iteration and 




Proof of convergence will follow naturally from the next discussion. Notice 
that this is a 2N x 2N system. We shall simplify it to an N x N system and 
increase the computational efficiency 

3.1 System Reduction 

We write L p and L™ in block matrix form as 



Mi -P 3 



-P 3 M 2 



pnp 4 _ F np 4 











where 



Mi = 2(aM{ + Mi) + Pi + P 3 
M 2 = 2M 2 + P 2 + P 3 

P = ^AtCj + At 2 K t 



It is easily seen that the left inverse of L p is given by 



-i 



D- 1 



-i 



where 



D = M l -P 3 M 2 1 P 3 

= 2aM f + [2Mi + Pi + P 3 (J - M^P Z )\ 

= {2aM f s + [2Mi + Pi + P 3 (J - M 2 ~ 1 P 3 )]W^- 1 }iy 



(3.15) 
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Note that D is invertible since Mf is positive definite, thus invertible, and 
all other terms are positive diagonal matrices, and thus their sum is positive 
definite and invertible. We also have 




-D- 1 T n P 4 



(3.16) 



Letting B n = (B[ l , B%), we have 



^n+l,fe+l _ 
-m+l,k+l _ 



W^iCi + D~ 1 Y n P i {u — v) n+1,k \ 
W- l M^Pz[Q + D^P^u - v) 



n+l,k 



(3.17) 
(3.18) 



where 



C^D^ + M^P^} 
<? = J D B - 1 [^ + MiP 3 - 1 SJ] 



(3.20) 



(3.19) 



At each time step, we do 2 iVx iV matrix solves in (3.19) and (3.20) to initialize 
the iterative scheme. Then, since the same term appears in both equations 
(3.17) and (3.18), for each k we only have to do 1 N x N matrix solve. In 
practice, since D s is symmetric, positive definite and time-independent, we 
compute the Cholesky factorization of D s at the start of the simulation and 
use the factorization for more efficient matrix solves at each step. As a side 
note, if we subtract (3.18) from (3.17), we have one equation for the OHC 
displacement u — v. 



4 Numerical Results 

4-1 Model Parameters 

We start with a modification of the parameters in [9] (See Table 1). It is known 
that higher dimensional models give higher sensitivity. This is the case with 
this model. The 1-D model [9] gives a 90 dB active gain at 16 kHz, whereas 
the 2-D model gives a 160 dB active gain. Thus, we need to tune the system to 
reduce the gain. There are many ways to do this, and the method we choose 
is to increase all the damping coefficients in the table by the following: 



2e' 



,0.27732 



i = l,2,3,4 
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Table 1 

Model parameters in cgs units 



m\{x) 


3 • 10~ 3 g • cur 2 


m m 


34.4 • 10 _ 


3 —2 

° g • cm 


ci(x) 


20 + 1500e" 2:r dyn • s • cm" 3 


Cm 


1.21 • 10 3 


dyn • 


s • cm 3 


fei(x) 


1.1 • 10 9 e" 4:c dyn-cm" 3 


h 


2.18 • 10 7 


dyn • 


cm~ 3 


m 2 (x) 


0.5 • 10~ 3 g • cm" 2 


L 


2.5 cm 






C2(x) 


10e~ 2 - 2x dyn • s • cm -3 


H 


0.1 cm 






k 2 (x) 


7 • 10 6 e~ 4Ax dyn • cm" 3 


P 


0.1 g • cm 


-3 




cs(x) 


2e~°- 8x dyn • s • cm~ 3 


e 


0.5 






k 3 (x) 


10 7 e~ 4x dyn • cm~ 3 


A 


0.08 cm 






c 4 (x) 


1040 e~ 2x dyn • s • cm~ 3 


At 


2.5 • 10~ 6 


- 10 


" 5 s 


k^{x) 


6.15 • W 8 e~ 4x dyn - cm" 3 


N 


401 
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Fig. 3. Both figures are sensitivity curves for CP = 0.77 cm or CF = 10 kHz. The 
left plot is a collection of sensitivity curves for the linear steady state active model 
where the parameter is the active gain 7. The right plot is a collection of sensitivity 
curves for the nonlinear time domain model where the parameter is pressure at the 
eardrum in dB SPL (sound pressure level). 



4-2 Isointensity Curves 



In an isocontour plot, a probe is placed at a specific location on the BM where 
the time response is measured and analyzed for input tones covering a range 
of frequencies. Figure 3 shows isointensity curves for CF = 10 kHz, which 
corresponds to CP = 0.77 cm. The characteristic place (CP) for a frequency 
is defined as the location on the BM of maximal response from a pure tone 
of that frequency in the fully linear active model (7 = 1). The characteristic 



10 




Time (ms) Time (ms) 

Fig. 4. An impulse, or click, lasting 0.1 ms starting at 0.4 ms is input into the 
nonlinear nonlocal ear model. The left plot is the BM displacement time series for 
various CF's ranging from 0.5-4 kHz. The right plot is a sensitivity plot for various 
stimulus intensities at CF = 6.4 kHz. 



frequency (CF) at a BM location is the inverse of this map. The left plot is 
the linear steady state active case. The parameter is the active gain 7, and 
for each value of the active gain we get a curve that is a function of the input 
frequency. The value of this function is the ratio |w|(CP)/P c , where |u|(CP) is 
BM displacement at the characteristic place and P e is pressure at the eardrum. 
This is known as sensitivity. It is basically an output /input ratio and gives 
the transfer characteristics of the ear at that particular active level. Notice 
that when 7 = 1, the BM at the characteristic place is most sensitive at the 
corresponding characteristic frequency, but at lower values of the gain, the 
sensitivity peak shifts to lower frequencies. 



Analogously, the second plot in Figure 3 shows isointensity curves for the non- 
linear time domain model where now the parameter is the intensity of the input 
stimulus in dB SPL (sound pressure level). For the time domain, we measure 
the root-mean-square BM amplitude from 5 ms (to remove transients) up to a 
certain time T. Note that for high-intensity tones, the model becomes passive 
while low-intensity tones give a more active model. This shows compression. 
Again, there is a frequency shift of the sensitivity peak (about one-half octave) 
from low to high-intensity stimuli in agreement with [11], so called half-octave 
shift. The plot agrees well with Figure 5 in [11]. 
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Time (ms) Time (ms) 

Fig. 5. Gaussian noise is input into the ear. The left plot is the BM displacement 
times series for various CF's ranging from 0.5-16 kHz. The right plot is a sensitivity 
plot for CF = 6.4 kHz. 

4-3 Complex Stimuli 

The first non-sinusoidal input we look at is a click. In the experiment in the 
left plot of Figure 4, we put probes at varying characteristic places associated 
with frequencies ranging from 0.5-4 kHz to measure the time series BM dis- 
placement. The click was 40 dB with duration 0.1 ms starting at 0.4 ms. All 
responses were normalized to amplitude 1. The plot is similar to Figure 4 in 
[5]. In the right plot of Figure 4, a probe was placed at CP for 6.4 kHz and the 
time series BM volume velocity was recorded for various intensities and the 
sensitivity plotted. This shows, similar to Figure 3, the compression effects at 
higher intensities. See Figure 9 in [11] for a similar plot. 

The second non-sinusoidal input we explore is Gaussian white noise. Figure 5 is 
similar in all regards to Figure 4. Notice again in the right plot the compression 
effect. 



4-4 Difference Tones 

Any nonlinear system with multiple sinusoidal inputs will create difference 
tones. If two frequencies f\ and f 2 are put into the ear, nfi±mf 2 will be created 
at varying intensities, where n and m are nonnegative integers. The cubic 
difference tone, denoted / = 2f\ — f 2 , where f\ < f 2 , is the most prominent. 
Figure 6 contains three plots of one experiment. The experiment consists of 
two sinusoidal tones, 7 and 10 kHz at 80 dB each. The cubic difference tone 
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BM location BM location Frequency (kHz) 

Fig. 6. Two sinusoidal tones, 7 and 10 kHz at 80 dB each, are the input. The 
left and middle plots are snapshots at 15 ms of BM displacement and active gain, 
respectively. The right plot is a spectrum plot of the BM displacement time series 
at CP for 4 kHz. 

is 4 kHz. The plot on the left is the BM profile for the experiment at 15 ms. 
We see combination tone peaks at 1.21 cm (CP for 4 kHz), 1.54 cm (CP for 
2 kHz) and 1.85 cm (CP for 1 kHz). The middle plot shows the snapshot at 
15 ms of the active gain parameter, showing the difference tones getting an 
active boost. Finally, the right plot is a spectrum plot of the time series for 
BM displacement at 1.21 cm, the characteristic place for 4 kHz. The cubic 
difference tone is above 1 nm and can therefore be heard. 



4-5 Multi-tone Suppression 

Two-tone (and multi-tone) suppression is characteristic of a compressive non- 
linearity and has been recognized in the ear [11,4,6]. Figure 7 illustrates two- 
tone suppression and is a collection of isodisplacement curves that show de- 
creased tuning in the presence of suppressors and is similar to Figure 16 in 
[11]. We placed a probe at the CP for 4 kHz (1.21 cm) and input sinusoids of 
various frequencies. At each frequency, we record the pressure at the eardrum 
that gives a 1 nm displacement for 4 kHz in the FFT spectrum of the time 
series response at CP. The curve without suppressors is dashed with circles. 
We then input each frequency again, but this time in the presence of a low 
side (0.5 kHz) tone and high side (7.5 kHz) tone, both at 80 dB. Notice the 
reduced tuning at the CF. Also notice the asymmetry of suppression, which 
shows low side is more suppressive than high side, in agreement with [6]. 

For multi-tone suppression, we look at tonal suppression of noise. In Figure 
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Frequency (kHz) 

Fig. 7. Iso displacement curves at CP for 4 kHz showing responses to single tones 
(dashed line w/ circle) and responses to the same tones in the presence of high-side 
and low-side suppressors presented at 80 dB SPL. 

8, for each plot, a probe was placed at every grid point along the BM and 
the time response was measured from 15 ms up to 25 ms. The signal in each 
consisted of noise at 50 dB with a 2 kHz tone ranging from 40 dB to 80 dB 
(top to bottom). An FFT was performed for each response and its charac- 
teristic frequency amplitude was recorded and plotted in decibels relative to 
the average of the response spectrum of dB noise from 0.5-16 kHz. We see 
suppression of all frequencies, with again low-side suppression stronger than 
high-side suppression. Figure 8 is qualitatively similar to Figure 3 in [4]. It is 
useful to compare this figure with Figure 9. This figure is the same as Figure 
8, except we do an FFT of the input signal at the eardrum. Comparing these 
two figures shows that we have a new spectral transform that can be used in 
place of an FFT in certain applications, for example signal recognition and 
noise suppression. 



5 Conclusions 



We studied a two-dimensional nonlinear nonlocal variation of the linear active 
model in [9]. We then developed an efficient and accurate numerical method 
and used this method to explore nonlinear effects of multi-tone sinusoidal 
inputs, as well as clicks and noise. We showed numerical results illustrating 
compression, multi-tone suppression and difference tones. The model reached 
agreement with experiments [11] and produced a novel nonlinear spectrum. 
In future work, we will analyze the model responses to speech and resulting 
spectra for speech recognition. It is also interesting to study the inverse prob- 
lem [12] of finding efficient and automated ways to tune the model to different 
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Fig. 8. Spectrum plots of BM responses for characteristic frequencies along the BM, 
from 500 Hz to 16 kHz, with 50 dB noise and a 2 kHz tone ranging from 40-80 dB. 
R0 is the average of the BM response spectrum of dB noise from 0.5-16 kHz. The 
solid line represents noise with tone, the dotted line noise without tone. 

physiological data. Applying the model to psychoacoustic signal processing 
[15] will be another fruitful line of inquiry. 
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A Appendix: Convergence of Iterative Scheme (3.14) 

We need the following Lemma: 
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Fig. 9. Spectrum plots of input signals consisting of 50 dB noise and a 2 kHz tone 
ranging from 40-80 dB. NO is the average of the spectrum of dB noise from 0.5-16 
kHz. The solid line represents noise with tone, the dotted line noise without tone. 

Lemma 1 If 

A -A 

M = 

B -B 

then every non-zero eigenvalue of M is an eigenvalue of A — B. 



PROOF. 

Let A be a non-zero eigenvalue of M with non-trivial eigenvector x = (xi,x 2 ). 
Thus, Mx = \x gives 

A(xi - x 2 ) = Xxi (A.l) 
B(x! - x 2 ) = \x 2 (A.2) 
Subtracting the two equations, we have 

(A - B)(xi - x 2 ) = A(fi - x 2 ) 

Now, if x\ — x 2 = 0, then from A.l and A.2 above and A 7^ 0, we have 
x\ = x 2 = 0. But this means x = 0, which is a contradiction. Thus, A is an 
eigenvalue of A — B with non-trivial eigenvector x\ — x 2 . □ 
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Theorem 2 There exists a constant C > such that if At < C , then 

P{L- 1 LD < 1 

where p is the spectral radius. Thus, the iterative scheme converges. 



PROOF. 

By the above lemma applied to (3.16), with constant 7, we have 

= 1 a[(I-M 2 1 P 3 )D- 1 P A ] 
where a denotes spectrum. Thus, we have 

piL-'Ll) < 7 ||(J - M^P^W-'D^P^ 

<i\\(i - m 2 1 p 3 )w-%\\d;%\\p 4 \\ 2 

Now, let (X,x) be the eigen-pair of D s with A the smallest eigenvalue and 
= I. Note that A > since D s is positive definite. Thus, we have 1/A is 
the largest eigenvalue of -D" 1 , which gives 

\\d;%<i/\ 

Thus, using the definition of D s from (3.15), we have 
A = x T D s x 

= x T {2aM* + [2Mi + P x + P 3 (I - M^P^W^jx 
>x T {[2M x + P 1 + P 3 (7 - M^P^W'^x 
> min{[2mi + p 1 +p 3 (l - m^ 1 ^)]™ -1 } 

where lowercase represents diagonal entries. The third line above follows from 
2aMf being positive definite. Finally, we have 

p(l p 1 l:)<7||(/-m 2 - 1 p 3 )^ 1 ||2||A7 1 ||2||P4|| 2 

max[(l - m2 1 p 3 )w~ 1 ] max(p 4 ) 

— 1 1 

min{[2mi + p 1 + p 3 (l - m 2 p 3 )\vj- x } 

For At small enough, we have convergence. □ 

With our parameters, for convergence it is sufficient that At < 0.0008. In 
practice, however, convergence is seen for At as large as 0.01. 
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